for (i in 1:length(params))
print(paste('Parameter:', names(params)[i], ' - Value:', params[[i]], '- Class:', class(params[[i]])))
## [1] "Parameter: Dataset - Value: CHD2_iPSCs_and_organoids_PublicRepo - Class: character"
## [1] "Parameter: SEFile - Value: /group/testa/Project/CHD2/BulkRNAseq/results/PublicRepo/2.DEA/day50/Output/Savings/day50CbO.SE_deseq2_HT.rds - Class: character"
## [1] "Parameter: DEAList - Value: /group/testa/Project/CHD2/BulkRNAseq/results/PublicRepo/2.DEA/day50/Output/Savings/day50CbO.DEAList_HT.rds - Class: character"
## [1] "Parameter: HT - Value: /group/testa/Project/CHD2/BulkRNAseq/results/PublicRepo/2.DEA/day50/Output/Savings/day50CbO.deseqHTvsWT.rds - Class: character"
## [1] "Parameter: SavingFolder - Value: /group/testa/Project/CHD2/BulkRNAseq/results/PublicRepo/6.Enrichments/day50/Output/Savings/ - Class: character"
## [1] "Parameter: FiguresFolder - Value: /group/testa/Project/CHD2/BulkRNAseq/results/PublicRepo/6.Enrichments/day50/Output/Figures/ - Class: character"
## [1] "Parameter: FDRthr - Value: 0.05 - Class: numeric"
## [1] "Parameter: logFCthr - Value: 0.55 - Class: numeric"
## [1] "Parameter: TopGO - Value: BP_MF_CC - Class: character"
## [1] "Parameter: GoEnTh - Value: 1 - Class: numeric"
## [1] "Parameter: GoPvalTh - Value: 0.05 - Class: numeric"
## [1] "Parameter: NbName - Value: TopGO_day50_HT - Class: character"
## [1] "Parameter: SaveImages - Value: FALSE - Class: logical"library(RNASeqBulkExploratory)
library(SummarizedExperiment)
library(tidyr)
library(AnnotationDbi)
library(org.Hs.eg.db)
library(topGO)
library(sechm)
library(ggplot2)
library(grid)
library(gridExtra)
library(RColorBrewer)
library(cowplot)
library(data.table)
source('/group/testa/Users/oliviero.leonardi/myProjects/CHD2/BulkRNAseq/ContainerHome/CHD2_organoids/NoGradientBarplots.R')Dataset <- params$Dataset
logFCthr <- params$logFCthr
FDRthr <- params$FDRthr
FdrTh <- FDRthr
logFcTh <- logFCthr
SavingFolder <- ifelse(is.null(params$SavingFolder), getwd(), params$SavingFolder)
FiguresFolder <- ifelse(is.null(params$FiguresFolder), getwd(), params$FiguresFolder)
if (dir.exists(SavingFolder) == FALSE) {
dir.create(SavingFolder, recursive=TRUE)
}#SE object coming from DEA, but not containing specific contrast results
SE_DEA <- readRDS(params$SEFile)
SE_DEA <- SE_DEA[rowData(SE_DEA)$GeneName != '', ]
rownames(SE_DEA) <- rowData(SE_DEA)$GeneName
# List with differential expression results (all time-points)
DEA <- readRDS(params$DEAList)colvector <- c("#5ec962", "#e95462", "#2c728e")
names(colvector) <- c('All', 'Up', 'Down')if(! identical(rownames(SE_DEA), row.names(DEA$HT$res))){
stop('Expression data in SE and results from differential espression analysis are inconsistent.')
}
SE_DEA <- mergeDeaSE(SE_DEA, DEA$HT$res, subsetRowDataCols=NULL,
logFcCol='log2FoldChange', FdrCol='padj') #specify
## Renaming " log2FoldChange " to "logFC"
## Renaming " padj " to "FDR"17673 genes in 14 samples have been testes for differential expression.
The following number of genes are identified as differentially expressed:
Imposing a threshold of 0.55 on the Log2FC and 0.05 on the FDR (as specified in parameters), 5601 genes are selected: 4905 up-regulated genes and 4894 down-regulated genes.
The results of the differential expression analysis are visualized by Volcano plot. An interactive version is included in the html (only genes with FDR < threshold), while a static version is saved.
plotVolcanoSE(SE=SE_DEA, FdrTh=FDRthr, logFcTh=logFCthr,
FdrCeil=1e-10, logFcCeil=4, Interactive = FALSE)## Warning: Removed 343 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 17531 rows containing missing values or values outside the scale range
## (`geom_text_repel()`).
## Warning: ggrepel: 100 unlabeled data points (too many overlaps). Consider increasing max.overlaps
plotVolcanoSE(SE=SE_DEA, FdrTh=FDRthr, logFcTh=logFCthr,
FdrCeil=1e-10, logFcCeil=4, Interactive = TRUE)Heatmaps for DEGs, showing scaled vst values.
DEGs <- dplyr::filter(data.frame(rowData(SE_DEA)), FDR < FDRthr & abs(logFC) > log2(logFCthr))
ScaledCols <- c('darkblue', "purple","white","lightgoldenrod1", 'goldenrod1')# sechm::sechm(SE_DEA, genes=DEGs$GeneName, assayName="vst", gaps_at="Genotype", show_rownames=FALSE,
# top_annotation=c('Genotype'), hmcols=ScaledCols, show_colnames=TRUE,
# do.scale=TRUE, breaks=0.85, column_title = "Scaled Vst Values")Gene ontology enrichment analysis is performed on the set of 5601 genes using TopGO with Fisher statistics and weight01 algorithm.
For each specified domain of the ontology:
I generate vectors for the gene universe, all modulated genes, up-regulated genes and down-regulated genes in the format required by TopGo.
GeneVectors <- topGOGeneVectors(SE_DEA, FdrTh=FDRthr, logFcTh=logFCthr)## Gene vector contains levels: 0,1
## Gene vector contains levels: 0,1
## Gene vector contains levels: 0,1
Therefore:
Then I set parameters according to the gene ontology domains to be evaluated. By default, Biological Process and Molecular Function domains are interrogated.
BpEval <- ifelse(length(grep('BP', params$TopGO))!=0, TRUE, FALSE)
MfEval <- ifelse(length(grep('MF', params$TopGO))!=0, TRUE, FALSE)
CcEval <- ifelse(length(grep('CC', params$TopGO))!=0, TRUE, FALSE)On the basis of the analysis settings, the enrichment for Biological Process IS performed.
# I generate a list that contains the association between each gene and the GO terms that are associated to it
BPannHT <- topGO::annFUN.org(whichOnto="BP", feasibleGenes=names(GeneVectors$DEGenes),
mapping="org.Hs.eg.db", ID="symbol") %>% inverseList()
# Wrapper function for topGO analysis (see helper file)
ResBPAllHT <- topGOResults(Genes=GeneVectors$DEGenes, gene2GO=BPannHT, ontology='BP',
desc=NULL, nodeSize=5, algorithm='weight01', statistic='fisher',
EnTh=params$GoEnTh, PvalTh=params$GoPvalTh, minTerms=15, geneTh=4,
saveRes=TRUE, fileName='BPAllHT', outDir=SavingFolder)
## Gene vector contains levels: 0,1
##
## Building most specific GOs .....
## ( 11741 GO terms found. )
##
## Build GO DAG topology ..........
## ( 15403 GO terms and 35047 relations. )
##
## Annotating nodes ...............
## ( 13971 genes annotated to the GO terms. )
##
## -- Weight01 Algorithm --
##
## the algorithm is scoring 7512 nontrivial nodes
## parameters:
## test statistic: fisher
##
## Level 19: 3 nodes to be scored (0 eliminated genes)
##
## Level 18: 8 nodes to be scored (0 eliminated genes)
##
## Level 17: 12 nodes to be scored (31 eliminated genes)
##
## Level 16: 29 nodes to be scored (54 eliminated genes)
##
## Level 15: 69 nodes to be scored (143 eliminated genes)
##
## Level 14: 122 nodes to be scored (357 eliminated genes)
##
## Level 13: 213 nodes to be scored (772 eliminated genes)
##
## Level 12: 370 nodes to be scored (1767 eliminated genes)
##
## Level 11: 663 nodes to be scored (3882 eliminated genes)
##
## Level 10: 953 nodes to be scored (6056 eliminated genes)
##
## Level 9: 1156 nodes to be scored (7727 eliminated genes)
##
## Level 8: 1143 nodes to be scored (9641 eliminated genes)
##
## Level 7: 1041 nodes to be scored (11150 eliminated genes)
##
## Level 6: 825 nodes to be scored (12270 eliminated genes)
##
## Level 5: 505 nodes to be scored (13030 eliminated genes)
##
## Level 4: 266 nodes to be scored (13499 eliminated genes)
##
## Level 3: 111 nodes to be scored (13667 eliminated genes)
##
## Level 2: 22 nodes to be scored (13759 eliminated genes)
##
## Level 1: 1 nodes to be scored (13808 eliminated genes)# Wrapper function for topGO analysis (see helper file)
ResBPDownHT <- topGOResults(Genes=GeneVectors$DEGenesDown, gene2GO=BPannHT, ontology='BP',
desc=NULL, nodeSize=5, algorithm='weight01', statistic='fisher',
EnTh=params$GoEnTh, PvalTh=params$GoPvalTh, minTerms=15, geneTh=4,
saveRes=TRUE, fileName='BPDownHT', outDir=paste0(SavingFolder))
## Gene vector contains levels: 0,1
##
## Building most specific GOs .....
## ( 11741 GO terms found. )
##
## Build GO DAG topology ..........
## ( 15403 GO terms and 35047 relations. )
##
## Annotating nodes ...............
## ( 13971 genes annotated to the GO terms. )
##
## -- Weight01 Algorithm --
##
## the algorithm is scoring 6238 nontrivial nodes
## parameters:
## test statistic: fisher
##
## Level 19: 2 nodes to be scored (0 eliminated genes)
##
## Level 18: 6 nodes to be scored (0 eliminated genes)
##
## Level 17: 9 nodes to be scored (26 eliminated genes)
##
## Level 16: 14 nodes to be scored (49 eliminated genes)
##
## Level 15: 45 nodes to be scored (120 eliminated genes)
##
## Level 14: 94 nodes to be scored (228 eliminated genes)
##
## Level 13: 158 nodes to be scored (651 eliminated genes)
##
## Level 12: 288 nodes to be scored (1586 eliminated genes)
##
## Level 11: 511 nodes to be scored (3689 eliminated genes)
##
## Level 10: 757 nodes to be scored (5865 eliminated genes)
##
## Level 9: 951 nodes to be scored (7463 eliminated genes)
##
## Level 8: 941 nodes to be scored (9375 eliminated genes)
##
## Level 7: 898 nodes to be scored (11035 eliminated genes)
##
## Level 6: 743 nodes to be scored (12201 eliminated genes)
##
## Level 5: 464 nodes to be scored (12990 eliminated genes)
##
## Level 4: 232 nodes to be scored (13493 eliminated genes)
##
## Level 3: 103 nodes to be scored (13666 eliminated genes)
##
## Level 2: 21 nodes to be scored (13757 eliminated genes)
##
## Level 1: 1 nodes to be scored (13808 eliminated genes)
# Selection on enrichment of at least 2 is implemented (also to avoid depleted categories). Then categories are ranked by PVal and all the ones with Pval < th are selected. If the number is < minTerms, othter terms are included to reach the minimum number. GOTable(ResBPDownHT$ResSel, maxGO=20)ResBPUpHT <- topGOResults(Genes=GeneVectors$DEGenesUp, gene2GO=BPannHT, ontology='BP',
desc=NULL, nodeSize=5, algorithm='weight01', statistic='fisher',
EnTh=params$GoEnTh, PvalTh=params$GoPvalTh, minTerms=15, geneTh=4,
saveRes=TRUE, fileName='BPUpHT', outDir=SavingFolder)
## Gene vector contains levels: 0,1
##
## Building most specific GOs .....
## ( 11741 GO terms found. )
##
## Build GO DAG topology ..........
## ( 15403 GO terms and 35047 relations. )
##
## Annotating nodes ...............
## ( 13971 genes annotated to the GO terms. )
##
## -- Weight01 Algorithm --
##
## the algorithm is scoring 6162 nontrivial nodes
## parameters:
## test statistic: fisher
##
## Level 19: 2 nodes to be scored (0 eliminated genes)
##
## Level 18: 3 nodes to be scored (0 eliminated genes)
##
## Level 17: 9 nodes to be scored (24 eliminated genes)
##
## Level 16: 23 nodes to be scored (41 eliminated genes)
##
## Level 15: 53 nodes to be scored (135 eliminated genes)
##
## Level 14: 93 nodes to be scored (318 eliminated genes)
##
## Level 13: 164 nodes to be scored (673 eliminated genes)
##
## Level 12: 278 nodes to be scored (1549 eliminated genes)
##
## Level 11: 493 nodes to be scored (3657 eliminated genes)
##
## Level 10: 721 nodes to be scored (5675 eliminated genes)
##
## Level 9: 914 nodes to be scored (7218 eliminated genes)
##
## Level 8: 935 nodes to be scored (9260 eliminated genes)
##
## Level 7: 900 nodes to be scored (10921 eliminated genes)
##
## Level 6: 730 nodes to be scored (12164 eliminated genes)
##
## Level 5: 468 nodes to be scored (12977 eliminated genes)
##
## Level 4: 248 nodes to be scored (13493 eliminated genes)
##
## Level 3: 105 nodes to be scored (13667 eliminated genes)
##
## Level 2: 22 nodes to be scored (13759 eliminated genes)
##
## Level 1: 1 nodes to be scored (13807 eliminated genes)
#dir.create(paste0(SavingFolder, 'TopGO/BPUp'), recursive=TRUE)
#GOAnnotation(ResBPUp$ResSel, GOdata=ResBPUp$GOdata, SavingFolder=paste0(SavingFolder, 'TopGO/BPUp'), keytype='SYMBOL')GOTable(ResBPUpHT$ResSel, maxGO=20)topGOBarplotAll(TopGOResAll=ResBPAllHT$ResSel, TopGOResDown=ResBPDownHT$ResSel, TopGOResUp = ResBPUpHT$ResSel,
terms=12, pvalTh=0.05, plotTitle=NULL, gradient = FALSE, cols = colvector)
## Warning in topGOBarplotAll(TopGOResAll = ResBPAllHT$ResSel, TopGOResDown = ResBPDownHT$ResSel, : Please make sure you specified names of the cols vector correctly.
## See `?topGOBarplotAll`On the basis of the analysis settings, the enrichment for Molecular Function IS performed.
# I generate a list that contains the association between each gene and the GO terms that are associated to it
MFannHT <- topGO::annFUN.org(whichOnto='MF', feasibleGenes=names(GeneVectors$DEGenes),
mapping='org.Hs.eg.db', ID='symbol') %>% inverseList()
# Wrapper function for topGO analysis (see helper file)
ResMFAllHT <- topGOResults(Genes=GeneVectors$DEGenes, gene2GO=MFannHT, ontology='MF',
desc=NULL, nodeSize=5, algorithm='weight01', statistic='fisher',
EnTh=params$GoEnTh, PvalTh=params$GoPvalTh, minTerms=15, geneTh=4,
saveRes=TRUE, fileName='MFAllHT', outDir=SavingFolder)
## Gene vector contains levels: 0,1
##
## Building most specific GOs .....
## ( 4161 GO terms found. )
##
## Build GO DAG topology ..........
## ( 4629 GO terms and 5968 relations. )
##
## Annotating nodes ...............
## ( 14383 genes annotated to the GO terms. )
##
## -- Weight01 Algorithm --
##
## the algorithm is scoring 1454 nontrivial nodes
## parameters:
## test statistic: fisher
##
## Level 12: 7 nodes to be scored (0 eliminated genes)
##
## Level 11: 19 nodes to be scored (0 eliminated genes)
##
## Level 10: 31 nodes to be scored (48 eliminated genes)
##
## Level 9: 75 nodes to be scored (247 eliminated genes)
##
## Level 8: 143 nodes to be scored (1458 eliminated genes)
##
## Level 7: 242 nodes to be scored (3639 eliminated genes)
##
## Level 6: 300 nodes to be scored (4571 eliminated genes)
##
## Level 5: 318 nodes to be scored (6463 eliminated genes)
##
## Level 4: 232 nodes to be scored (9395 eliminated genes)
##
## Level 3: 69 nodes to be scored (11655 eliminated genes)
##
## Level 2: 17 nodes to be scored (12371 eliminated genes)
##
## Level 1: 1 nodes to be scored (14249 eliminated genes)ResMFDownHT <- topGOResults(Genes=GeneVectors$DEGenesDown, gene2GO=MFannHT, ontology='MF',
desc=NULL, nodeSize=5, algorithm='weight01', statistic='fisher',
EnTh=params$GoEnTh, PvalTh=params$GoPvalTh, minTerms=15, geneTh=4,
saveRes=TRUE, fileName='MFDownHT', outDir=SavingFolder)
## Gene vector contains levels: 0,1
##
## Building most specific GOs .....
## ( 4161 GO terms found. )
##
## Build GO DAG topology ..........
## ( 4629 GO terms and 5968 relations. )
##
## Annotating nodes ...............
## ( 14383 genes annotated to the GO terms. )
##
## -- Weight01 Algorithm --
##
## the algorithm is scoring 1184 nontrivial nodes
## parameters:
## test statistic: fisher
##
## Level 12: 3 nodes to be scored (0 eliminated genes)
##
## Level 11: 13 nodes to be scored (0 eliminated genes)
##
## Level 10: 19 nodes to be scored (29 eliminated genes)
##
## Level 9: 61 nodes to be scored (200 eliminated genes)
##
## Level 8: 112 nodes to be scored (1351 eliminated genes)
##
## Level 7: 189 nodes to be scored (3574 eliminated genes)
##
## Level 6: 234 nodes to be scored (4419 eliminated genes)
##
## Level 5: 264 nodes to be scored (6210 eliminated genes)
##
## Level 4: 208 nodes to be scored (9123 eliminated genes)
##
## Level 3: 63 nodes to be scored (11536 eliminated genes)
##
## Level 2: 17 nodes to be scored (12355 eliminated genes)
##
## Level 1: 1 nodes to be scored (14246 eliminated genes)
#dir.create(paste0(SavingFolder, 'TopGO/MFDown'), recursive=TRUE)
#GOAnnotation(ResMFDown$ResSel, GOdata=ResMFDown$GOdata, SavingFolder=paste0(SavingFolder, 'TopGO/MFDown'), keytype='SYMBOL')GOTable(ResMFDownHT$ResSel, maxGO=20)ResMFUpHT <- topGOResults(Genes=GeneVectors$DEGenesUp, gene2GO=MFannHT, ontology='MF',
desc=NULL, nodeSize=5, algorithm='weight01', statistic='fisher',
EnTh=params$GoEnTh, PvalTh=params$GoPvalTh, minTerms=15, geneTh=4,
saveRes=TRUE, fileName='MFUpHT', outDir=SavingFolder)
## Gene vector contains levels: 0,1
##
## Building most specific GOs .....
## ( 4161 GO terms found. )
##
## Build GO DAG topology ..........
## ( 4629 GO terms and 5968 relations. )
##
## Annotating nodes ...............
## ( 14383 genes annotated to the GO terms. )
##
## -- Weight01 Algorithm --
##
## the algorithm is scoring 1130 nontrivial nodes
## parameters:
## test statistic: fisher
##
## Level 12: 5 nodes to be scored (0 eliminated genes)
##
## Level 11: 17 nodes to be scored (0 eliminated genes)
##
## Level 10: 24 nodes to be scored (29 eliminated genes)
##
## Level 9: 52 nodes to be scored (236 eliminated genes)
##
## Level 8: 106 nodes to be scored (1388 eliminated genes)
##
## Level 7: 171 nodes to be scored (3537 eliminated genes)
##
## Level 6: 227 nodes to be scored (4397 eliminated genes)
##
## Level 5: 252 nodes to be scored (6115 eliminated genes)
##
## Level 4: 198 nodes to be scored (9100 eliminated genes)
##
## Level 3: 61 nodes to be scored (11574 eliminated genes)
##
## Level 2: 16 nodes to be scored (12342 eliminated genes)
##
## Level 1: 1 nodes to be scored (14231 eliminated genes)
#dir.create(paste0(SavingFolder, 'TopGO/MFUp'), recursive=TRUE)
#GOAnnotation(ResMFUp$ResSel, GOdata=ResMFUp$GOdata, SavingFolder=paste0(SavingFolder, 'TopGO/MFUp'), keytype='SYMBOL')GOTable(ResMFUpHT$ResSel, maxGO=20)topGOBarplotAll(TopGOResAll=ResMFAllHT$ResSel, TopGOResDown=ResMFDownHT$ResSel, TopGOResUp=ResMFUpHT$ResSel,
terms=12, pvalTh=0.05, plotTitle=NULL, gradient = FALSE, cols = colvector)
## Warning in topGOBarplotAll(TopGOResAll = ResMFAllHT$ResSel, TopGOResDown = ResMFDownHT$ResSel, : Please make sure you specified names of the cols vector correctly.
## See `?topGOBarplotAll`On the basis of the analysis settings, the enrichment for Cellular Component IS performed.
# I generate a list that contains the association between each gene and the GO terms that are associated to it
CCannHT <- topGO::annFUN.org(whichOnto='CC', feasibleGenes=names(GeneVectors$DEGenes),
mapping='org.Hs.eg.db', ID='symbol') %>% inverseList()
# Wrapper function for topGO analysis (see helper file)
ResCCAllHT <- topGOResults(Genes=GeneVectors$DEGenes, gene2GO=CCannHT, ontology='CC',
desc=NULL, nodeSize=5, algorithm='weight01', statistic='fisher',
EnTh=params$GoEnTh, PvalTh=params$GoPvalTh, minTerms=15, geneTh=4,
saveRes=TRUE, fileName='CCAllHT', outDir=SavingFolder)
## Gene vector contains levels: 0,1
##
## Building most specific GOs .....
## ( 1698 GO terms found. )
##
## Build GO DAG topology ..........
## ( 1919 GO terms and 3235 relations. )
##
## Annotating nodes ...............
## ( 14645 genes annotated to the GO terms. )
##
## -- Weight01 Algorithm --
##
## the algorithm is scoring 818 nontrivial nodes
## parameters:
## test statistic: fisher
##
## Level 13: 2 nodes to be scored (0 eliminated genes)
##
## Level 12: 6 nodes to be scored (0 eliminated genes)
##
## Level 11: 35 nodes to be scored (47 eliminated genes)
##
## Level 10: 86 nodes to be scored (77 eliminated genes)
##
## Level 9: 125 nodes to be scored (822 eliminated genes)
##
## Level 8: 130 nodes to be scored (2908 eliminated genes)
##
## Level 7: 129 nodes to be scored (5313 eliminated genes)
##
## Level 6: 113 nodes to be scored (8895 eliminated genes)
##
## Level 5: 81 nodes to be scored (10622 eliminated genes)
##
## Level 4: 53 nodes to be scored (12694 eliminated genes)
##
## Level 3: 55 nodes to be scored (13996 eliminated genes)
##
## Level 2: 2 nodes to be scored (14423 eliminated genes)
##
## Level 1: 1 nodes to be scored (14570 eliminated genes)
#write.table(ResCCAll$ResAll, file=paste0(SavingFolder, 'TopGO/CCAllResults.txt'), sep='\t', row.names=FALSE)# Wrapper function for topGO analysis (see helper file)
ResCCDownHT <- topGOResults(Genes=GeneVectors$DEGenesDown, gene2GO=CCannHT, ontology='CC',
desc=NULL, nodeSize=5, algorithm='weight01', statistic='fisher',
EnTh=params$GoEnTh, PvalTh=params$GoPvalTh, minTerms=15, geneTh=4,
saveRes=TRUE, fileName='CCDownHT', outDir=SavingFolder)
## Gene vector contains levels: 0,1
##
## Building most specific GOs .....
## ( 1698 GO terms found. )
##
## Build GO DAG topology ..........
## ( 1919 GO terms and 3235 relations. )
##
## Annotating nodes ...............
## ( 14645 genes annotated to the GO terms. )
##
## -- Weight01 Algorithm --
##
## the algorithm is scoring 676 nontrivial nodes
## parameters:
## test statistic: fisher
##
## Level 13: 1 nodes to be scored (0 eliminated genes)
##
## Level 12: 3 nodes to be scored (0 eliminated genes)
##
## Level 11: 26 nodes to be scored (35 eliminated genes)
##
## Level 10: 70 nodes to be scored (65 eliminated genes)
##
## Level 9: 107 nodes to be scored (761 eliminated genes)
##
## Level 8: 110 nodes to be scored (2822 eliminated genes)
##
## Level 7: 103 nodes to be scored (5220 eliminated genes)
##
## Level 6: 96 nodes to be scored (8853 eliminated genes)
##
## Level 5: 66 nodes to be scored (10575 eliminated genes)
##
## Level 4: 42 nodes to be scored (12676 eliminated genes)
##
## Level 3: 49 nodes to be scored (13994 eliminated genes)
##
## Level 2: 2 nodes to be scored (14422 eliminated genes)
##
## Level 1: 1 nodes to be scored (14570 eliminated genes)
#dir.create(paste0(SavingFolder, 'TopGO/CCDown'), recursive=TRUE)
#GOAnnotation(ResCCDown$ResSel, GOdata=ResCCDown$GOdata, SavingFolder=paste0(SavingFolder, 'TopGO/CCDown'), keytype='SYMBOL')GOTable(ResCCDownHT$ResSel, maxGO=20)# Wrapper function for topGO analysis (see helper file)
ResCCUpHT <- topGOResults(Genes=GeneVectors$DEGenesUp, gene2GO=CCannHT, ontology='CC',
desc=NULL, nodeSize=5, algorithm='weight01', statistic='fisher',
EnTh=params$GoEnTh, PvalTh=params$GoPvalTh, minTerms=15, geneTh=4,
saveRes=TRUE, fileName='CCUpHT', outDir=SavingFolder)
## Gene vector contains levels: 0,1
##
## Building most specific GOs .....
## ( 1698 GO terms found. )
##
## Build GO DAG topology ..........
## ( 1919 GO terms and 3235 relations. )
##
## Annotating nodes ...............
## ( 14645 genes annotated to the GO terms. )
##
## -- Weight01 Algorithm --
##
## the algorithm is scoring 664 nontrivial nodes
## parameters:
## test statistic: fisher
##
## Level 13: 2 nodes to be scored (0 eliminated genes)
##
## Level 12: 5 nodes to be scored (0 eliminated genes)
##
## Level 11: 26 nodes to be scored (47 eliminated genes)
##
## Level 10: 66 nodes to be scored (69 eliminated genes)
##
## Level 9: 92 nodes to be scored (761 eliminated genes)
##
## Level 8: 100 nodes to be scored (2623 eliminated genes)
##
## Level 7: 105 nodes to be scored (5044 eliminated genes)
##
## Level 6: 95 nodes to be scored (8743 eliminated genes)
##
## Level 5: 76 nodes to be scored (10591 eliminated genes)
##
## Level 4: 49 nodes to be scored (12677 eliminated genes)
##
## Level 3: 45 nodes to be scored (13996 eliminated genes)
##
## Level 2: 2 nodes to be scored (14422 eliminated genes)
##
## Level 1: 1 nodes to be scored (14570 eliminated genes)
#dir.create(paste0(SavingFolder, 'TopGO/CCUp'), recursive=TRUE)
#GOAnnotation(ResCCUp$ResSel, GOdata=ResCCUp$GOdata, SavingFolder=paste0(SavingFolder, 'TopGO/CCUp'), keytype='SYMBOL')GOTable(ResCCUpHT$ResSel, maxGO=20)topGOBarplotAll(TopGOResAll=ResCCAllHT$ResSel, TopGOResDown=ResCCDownHT$ResSel, TopGOResUp=ResCCUpHT$ResSel,
terms=12, pvalTh=0.05, plotTitle=NULL, gradient = FALSE, cols = colvector)
## Warning in topGOBarplotAll(TopGOResAll = ResCCAllHT$ResSel, TopGOResDown = ResCCDownHT$ResSel, : Please make sure you specified names of the cols vector correctly.
## See `?topGOBarplotAll`Most of the useful information has been saved during the analysis. Here I save figures, workspace and information about the session.
if (params$SaveImages == TRUE){ #Just in case since eval only works when knitting
#Set the folder paths
from <- paste(getwd(), paste(params$NbName, 'files/figure-html', sep='_'), sep='/')
to <- params$FiguresFolder
#Copy to output directory
file.copy(from, to, recursive = TRUE, copy.mode = TRUE)
}ResSelBP_HT <- list(ResSelAll = ResBPAllHT$ResSel,
ResSelUp = ResBPUpHT$ResSel,
ResSelDown = ResBPDownHT$ResSel)
ResSelMF_HT <- list(ResSelAll = ResMFAllHT$ResSel,
ResSelUp = ResMFUpHT$ResSel,
ResSelDown = ResMFDownHT$ResSel)
ResSelCC_HT <- list(ResSelAll = ResCCAllHT$ResSel,
ResSelUp = ResCCUpHT$ResSel,
ResSelDown = ResCCDownHT$ResSel)
ResSel_HT <- list(BP = ResSelBP_HT,
MF = ResSelMF_HT,
CC = ResSelCC_HT)
saveRDS(ResSel_HT, paste0(SavingFolder, '/day50CbO.', 'ResSel_HT.rds'))SessionInfo <- sessionInfo()
Date <- date()
save.image(paste0(SavingFolder, '/day50CbO.', 'FunctionalAnalysisWorkspace_HT.RData'))